function [residual, g1, g2, g3] = a15_markups_dynamic(y, x, params, steady_state, it_)
%
% Status : Computes dynamic model for Dynare
%
% Inputs :
%   y         [#dynamic variables by 1] double    vector of endogenous variables in the order stored
%                                                 in M_.lead_lag_incidence; see the Manual
%   x         [M_.exo_nbr by nperiods] double     matrix of exogenous variables (in declaration order)
%                                                 for all simulation periods
%   params    [M_.param_nbr by 1] double          vector of parameter values in declaration order
%   it_       scalar double                       time period for exogenous variables for which to evaluate the model
%
% Outputs:
%   residual  [M_.endo_nbr by 1] double    vector of residuals of the dynamic model equations in order of 
%                                          declaration of the equations
%   g1        [M_.endo_nbr by #dynamic variables] double    Jacobian matrix of the dynamic model equations;
%                                                           rows: equations in order of declaration
%                                                           columns: variables in order stored in M_.lead_lag_incidence
%   g2        [M_.endo_nbr by (#dynamic variables)^2] double   Hessian matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence
%   g3        [M_.endo_nbr by (#dynamic variables)^3] double   Third order derivative matrix of the dynamic model equations;
%                                                              rows: equations in order of declaration
%                                                              columns: variables in order stored in M_.lead_lag_incidence
%
%
% Warning : this file is generated automatically by Dynare
%           from model file (.mod)

%
% Model equations
%

residual = zeros(21, 1);
T47 = (exp(y(7))^(1-params(11))+params(10)-1)/params(10);
T51 = T47^((params(11)-params(7))/(1-params(11)));
T99 = params(1)*exp(y(33))*exp(y(31)-y(13))^2*y(35);
T150 = 1+(params(10)-1)*exp(y(7))^(params(11)-1);
T230 = (exp(y(7))*exp(y(19))*(exp(y(7))*exp(y(19))-1)-exp(y(7))*exp(y(19))*exp(y(7))*exp(y(19)))/((exp(y(7))*exp(y(19))-1)*(exp(y(7))*exp(y(19))-1));
lhs =y(21);
rhs =params(13)/2*(exp(y(13)-y(3))-1)^2;
residual(1)= lhs-rhs;
lhs =y(22);
rhs =params(13)*(exp(y(13)-y(3))-1);
residual(2)= lhs-rhs;
lhs =y(24);
rhs =y(29)-y(10);
residual(3)= lhs-rhs;
lhs =y(20);
rhs =log(params(1))+y(32)-y(15);
residual(4)= lhs-rhs;
lhs =exp(y(8));
rhs =exp(y(7))^(-params(11))*T51;
residual(5)= lhs-rhs;
lhs =exp(y(15));
rhs =(exp(y(11))-params(14)*exp(y(1)))^(-params(8))-params(1)*params(14)*(exp(y(30))-exp(y(11))*params(14))^(-params(8));
residual(6)= lhs-rhs;
lhs =params(9)*exp(y(14))^params(2);
rhs =exp(y(15))*exp(y(18));
residual(7)= lhs-rhs;
lhs =exp(y(15));
rhs =exp(y(16))*(1-y(21)-exp(y(13)-y(3))*y(22))+T99;
residual(8)= lhs-rhs;
lhs =exp(y(16));
rhs =params(1)*(exp(y(32))*exp(y(34))+exp(y(33))*(1-params(3)));
residual(9)= lhs-rhs;
lhs =(1-y(21))*exp(y(13));
rhs =exp(y(12))-(1-params(3))*exp(y(2));
residual(10)= lhs-rhs;
lhs =y(10);
rhs =y(25)+y(2)*params(4)+y(14)*(1-params(4));
residual(11)= lhs-rhs;
lhs =y(17);
rhs =y(10)+(-y(19))+log(params(4))-y(2);
residual(12)= lhs-rhs;
lhs =y(18);
rhs =y(10)+(-y(19))+log(1-params(4))-y(14);
residual(13)= lhs-rhs;
lhs =exp(y(7))*exp(y(19))/(exp(y(7))*exp(y(19))-1);
rhs =params(11)-(params(11)-params(7))/T150;
residual(14)= lhs-rhs;
lhs =exp(y(8))*(exp(y(7))-exp((-y(19))))-(1-exp((-y(19))));
rhs =params(5)*exp(y(20))*exp(y(24)+y(20)+y(28));
residual(15)= lhs-rhs;
lhs =exp(y(9));
rhs =1-exp((-y(19)))-(1-1/params(12))*(params(12)/exp(y(19)))^(1-params(7))+params(5)*exp(y(24)+y(20)+y(28));
residual(16)= lhs-rhs;
lhs =exp(y(11))+exp(y(13))+exp(y(23));
rhs =exp(y(10));
residual(17)= lhs-rhs;
lhs =y(23);
rhs =params(6)+y(26);
residual(18)= lhs-rhs;
lhs =y(25);
rhs =params(15)*y(4)+x(it_, 1);
residual(19)= lhs-rhs;
lhs =y(26);
rhs =x(it_, 2)+params(16)*y(5)+params(17)*y(6);
residual(20)= lhs-rhs;
lhs =y(27);
rhs =y(5);
residual(21)= lhs-rhs;
if nargout >= 2,
  g1 = zeros(21, 37);

  %
  % Jacobian matrix
  %

  g1(1,3)=(-(params(13)/2*(-exp(y(13)-y(3)))*2*(exp(y(13)-y(3))-1)));
  g1(1,13)=(-(params(13)/2*exp(y(13)-y(3))*2*(exp(y(13)-y(3))-1)));
  g1(1,21)=1;
  g1(2,3)=(-(params(13)*(-exp(y(13)-y(3)))));
  g1(2,13)=(-(params(13)*exp(y(13)-y(3))));
  g1(2,22)=1;
  g1(3,10)=1;
  g1(3,29)=(-1);
  g1(3,24)=1;
  g1(4,15)=1;
  g1(4,32)=(-1);
  g1(4,20)=1;
  g1(5,7)=(-(T51*exp(y(7))*getPowerDeriv(exp(y(7)),(-params(11)),1)+exp(y(7))^(-params(11))*exp(y(7))*getPowerDeriv(exp(y(7)),1-params(11),1)/params(10)*getPowerDeriv(T47,(params(11)-params(7))/(1-params(11)),1)));
  g1(5,8)=exp(y(8));
  g1(6,1)=(-((-(params(14)*exp(y(1))))*getPowerDeriv(exp(y(11))-params(14)*exp(y(1)),(-params(8)),1)));
  g1(6,11)=(-(exp(y(11))*getPowerDeriv(exp(y(11))-params(14)*exp(y(1)),(-params(8)),1)-params(1)*params(14)*(-(exp(y(11))*params(14)))*getPowerDeriv(exp(y(30))-exp(y(11))*params(14),(-params(8)),1)));
  g1(6,30)=params(1)*params(14)*exp(y(30))*getPowerDeriv(exp(y(30))-exp(y(11))*params(14),(-params(8)),1);
  g1(6,15)=exp(y(15));
  g1(7,14)=params(9)*exp(y(14))*getPowerDeriv(exp(y(14)),params(2),1);
  g1(7,15)=(-(exp(y(15))*exp(y(18))));
  g1(7,18)=(-(exp(y(15))*exp(y(18))));
  g1(8,3)=(-(exp(y(16))*(-(y(22)*(-exp(y(13)-y(3)))))));
  g1(8,13)=(-(exp(y(16))*(-(exp(y(13)-y(3))*y(22)))+y(35)*params(1)*exp(y(33))*(-exp(y(31)-y(13)))*2*exp(y(31)-y(13))));
  g1(8,31)=(-(y(35)*params(1)*exp(y(33))*exp(y(31)-y(13))*2*exp(y(31)-y(13))));
  g1(8,15)=exp(y(15));
  g1(8,16)=(-(exp(y(16))*(1-y(21)-exp(y(13)-y(3))*y(22))));
  g1(8,33)=(-T99);
  g1(8,21)=exp(y(16));
  g1(8,22)=(-(exp(y(16))*(-exp(y(13)-y(3)))));
  g1(8,35)=(-(params(1)*exp(y(33))*exp(y(31)-y(13))^2));
  g1(9,32)=(-(params(1)*exp(y(32))*exp(y(34))));
  g1(9,16)=exp(y(16));
  g1(9,33)=(-(params(1)*exp(y(33))*(1-params(3))));
  g1(9,34)=(-(params(1)*exp(y(32))*exp(y(34))));
  g1(10,2)=(1-params(3))*exp(y(2));
  g1(10,12)=(-exp(y(12)));
  g1(10,13)=(1-y(21))*exp(y(13));
  g1(10,21)=(-exp(y(13)));
  g1(11,10)=1;
  g1(11,2)=(-params(4));
  g1(11,14)=(-(1-params(4)));
  g1(11,25)=(-1);
  g1(12,10)=(-1);
  g1(12,2)=1;
  g1(12,17)=1;
  g1(12,19)=1;
  g1(13,10)=(-1);
  g1(13,14)=1;
  g1(13,18)=1;
  g1(13,19)=1;
  g1(14,7)=T230-(-((-((params(11)-params(7))*(params(10)-1)*exp(y(7))*getPowerDeriv(exp(y(7)),params(11)-1,1)))/(T150*T150)));
  g1(14,19)=T230;
  g1(15,7)=exp(y(8))*exp(y(7));
  g1(15,8)=exp(y(8))*(exp(y(7))-exp((-y(19))));
  g1(15,28)=(-(params(5)*exp(y(20))*exp(y(24)+y(20)+y(28))));
  g1(15,19)=exp(y(8))*exp((-y(19)))-exp((-y(19)));
  g1(15,20)=(-(params(5)*exp(y(20))*exp(y(24)+y(20)+y(28))+params(5)*exp(y(20))*exp(y(24)+y(20)+y(28))));
  g1(15,24)=(-(params(5)*exp(y(20))*exp(y(24)+y(20)+y(28))));
  g1(16,9)=exp(y(9));
  g1(16,28)=(-(params(5)*exp(y(24)+y(20)+y(28))));
  g1(16,19)=(-(exp((-y(19)))-(1-1/params(12))*(-(exp(y(19))*params(12)))/(exp(y(19))*exp(y(19)))*getPowerDeriv(params(12)/exp(y(19)),1-params(7),1)));
  g1(16,20)=(-(params(5)*exp(y(24)+y(20)+y(28))));
  g1(16,24)=(-(params(5)*exp(y(24)+y(20)+y(28))));
  g1(17,10)=(-exp(y(10)));
  g1(17,11)=exp(y(11));
  g1(17,13)=exp(y(13));
  g1(17,23)=exp(y(23));
  g1(18,23)=1;
  g1(18,26)=(-1);
  g1(19,4)=(-params(15));
  g1(19,25)=1;
  g1(19,36)=(-1);
  g1(20,5)=(-params(16));
  g1(20,26)=1;
  g1(20,37)=(-1);
  g1(20,6)=(-params(17));
  g1(21,5)=(-1);
  g1(21,27)=1;
end
if nargout >= 3,
  %
  % Hessian matrix
  %

  g2 = sparse([],[],[],21,1369);
end
if nargout >= 4,
  %
  % Third order derivatives
  %

  g3 = sparse([],[],[],21,50653);
end
end
